﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace CM1
{
    class Spline
    {   // Коэффициенты сплайна
        // Канонический вид сплайна 
        // a_i + b_i * (x - x_i) + c_i * (x - x_i)^2 / 2.0 + d_i * (x - x_i)^3 / 6.0
        // i = 1 .. n - 1
        public double[] a;
        public double[] b;
        public double[] c;
        public double[] d;

        // Число узлов сетки
        public int n;

        //Значения производных сплайна на границе
        public double LP, RP;

        // Значения аргументов для интерполяции
        public double[] x;

        // Значения функции в узлах сетки
        public double[] y;

        // Вычисление значения аппроксимируемой функции и её первой производной
        // функция
        public double ComputeFunc(double arg, int type)
        {
            switch (type)
            {
                case 1:
                    if (arg <= 0.0) return arg * arg * arg + 3 * arg * arg;
                    else return -arg * arg * arg + 3 * arg * arg;
                case 2:
                    return System.Math.Sqrt(arg) * System.Math.Cos(arg);
                case 3:
                    return System.Math.Sqrt(arg) * System.Math.Cos(arg) + System.Math.Cos(10 * arg);
                case 4:
                    return System.Math.Sqrt(arg) * System.Math.Cos(arg) + System.Math.Cos(100 * arg);
                default: return 999999999999;
            }
        }
        // производная
        public double ComputePrFunc(double arg, int type)
        {
            switch (type)
            {
                case 1:
                    if (arg <= 0.0) return 3 * arg * arg + 6 * arg;
                    else return -3 * arg * arg + 6 * arg;
                case 2:
                    return System.Math.Cos(arg) / (2 * System.Math.Sqrt(arg))
                        - System.Math.Sin(arg) * System.Math.Sqrt(arg);
                case 3:
                    return System.Math.Cos(arg) / (2.0 * System.Math.Sqrt(arg))
                        - System.Math.Sin(arg) * System.Math.Sqrt(arg) - 10.0 * System.Math.Sin(10 * arg);
                case 4:
                    return System.Math.Cos(arg) / (2 * System.Math.Sqrt(arg))
                        - System.Math.Sin(arg) * System.Math.Sqrt(arg) - 100.0 * System.Math.Sin(100 * arg);

                default: return 999999999999;
            }
        }
        // Вставить в сплайн исходные данные
        public void CSpline(int n1, int GUtype, int type)
        {
            // Число узлов
            n = n1;

            // Выделение памяти под сплайн
            a = new double[n];
            b = new double[n];
            c = new double[n];
            d = new double[n];
            x = new double[n];
            y = new double[n];

            // Вычислим ГУ
            if (GUtype == 1)
            {
                // ЕГУ
                LP = 0.0;
                RP = 0.0;
            }
            if (GUtype == 2)
            {
                if (type == 1)
                {
                    LP = 0.0;
                    RP = 0.0;
                }
                if (type == 2)
                {
                    LP = -1.5168488671434721;
                    RP = 1.817350631436782;
                }
                if (type == 3)
                {
                    LP = 82.39030404049574;
                    RP = -98.18264936856626;
                }
                if (type == 4)
                {
                    LP = -8624.705571858944;
                    RP = -9998.182649174905;
                }
            }

            // Вычисление аргументов узлов и значений фукнции в соответсвующих точках
            if (type == 1)
            {
                for (int i = 0; i < n; i++)
                {

                    x[i] = -1.0 + i * 2.0 / (n - 1.0);
                    y[i] = ComputeFunc(x[i], 1);
                }
            }
            else
            {
                for (int i = 0; i < n; i++)
                {

                    x[i] = 1.0 + i * 2.1415926535897932384626433832795 / (n - 1.0);
                    y[i] = ComputeFunc(x[i], type);
                }
            }
        }
        public void ComputeSpline()
        {
            // Инициализация массива сплайнов
            for (int i = 0; i < n; ++i)
            {
                a[i] = y[i];
            }
            c[0] = LP;
            c[n - 1] = RP;

            // Вычисление коэффициентов c[i] методом прогонки
            // Прямой Ход
            double[] alpha = new double[n - 1];
            double[] beta = new double[n - 1];
            alpha[0] = 0;
            beta[0] = LP;
            for (int i = 1; i < n - 1; ++i)
            {
                double h_i = x[i] - x[i - 1], h_i1 = x[i + 1] - x[i];
                double A = h_i;
                double C = 2.0 * (h_i + h_i1);
                double B = h_i1;
                double F = 6.0 * ((y[i + 1] - y[i]) / h_i1 - (y[i] - y[i - 1]) / h_i);
                double z = (A * alpha[i - 1] + C);
                alpha[i] = -B / z;
                beta[i] = (F - A * beta[i - 1]) / z;
            }

            //Обратный ход прогонки - вычисление коэфф c_i
            for (int i = n - 2; i > 0; --i)
                c[i] = alpha[i] * c[i + 1] + beta[i];

            // Очистка памяти
            alpha = null;
            beta = null;

            //Далее находим a[i] b[i] d[i]
            for (int i = n - 1; i > 0; --i)
            {
                double h_i = x[i] - x[i - 1];
                d[i] = (c[i] - c[i - 1]) / h_i;
                b[i] = h_i * c[i] / 2.0 - h_i * h_i * d[i] / (6.0) + (y[i] - y[i - 1]) / h_i;
            }
        }
        public double CompSplineVal(double arg)
        {
            int s = -1;
            if (arg <= x[0]) // Если x меньше точки сетки x[0] - пользуемся первым сплайном
                s = 1;
            else if (arg >= x[n - 1]) // Если x больше точки сетки x[n - 1] - пользуемся последним сплайном
                s = n - 1;
            else // Иначе x лежит между граничными точками сетки - производим поиск нужного сплайна
            {
                for (int i = 1; i < n; i++)
                {
                    if (arg < x[i])
                    {
                        s = i;
                        break;
                    }
                }
            }

            double dx = (arg - x[s]);
            return a[s] + b[s] * dx + c[s] * dx * dx / 2.0 + d[s] * dx * dx * dx / 6.0;
        }
        // Вычисление значения производной сплайна в точке
        public double CompSplineValPr(double arg)
        {
            int s = -1;
            if (arg <= x[0]) // Если x меньше точки сетки x[0] - пользуемся первым сплайном
                s = 1;
            else if (arg >= x[n - 1]) // Если x больше точки сетки x[n - 2] - пользуемся последним сплайном
                s = n - 1;
            else // Иначе x лежит между граничными точками сетки - производим поиск нужного сплайна
            {
                for (int i = 1; i < n; i++)
                {
                    if (arg < x[i])
                    {
                        s = i;
                        break;
                    }
                }
            }

            double dx = (arg - x[s]);
            return b[s] + c[s] * dx + d[s] * dx * dx / 2.0;
        }
    }
}
